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V.I. Usachenko and S.-I. Chu [Phys. Rev. A 71, 063410 (2005)] discuss the molecular strong-field 
approximation in the velocity gauge formulation and indicate that some of our earlier velocity gauge 
calculations are inaccurate. Here we comment on the results of Usachenko and Chu. First, we show 
that the molecular orbitals used by Usachenko and Chu do not have the correct symmetry, and 
second, that it is an oversimplification to describe the molecular orbitals in terms of just a single 
linear combination of two atomic orbitals. Finally, some values for the generalized Bessel function 
are given for comparison. 

PACS numbers: 33.80.Rv,32.80.Rm 



I. INTRODUCTION 



II. MOLECULAR ORBITALS 



The molecular strong-field approximation (MO-SFA) 
was formulated in the velocity gauge in Ref. [1], and more 
recently in the length gauge Q|. The latter work shows 
that the velocity gauge MO-SFA with the initial state of 
the highest occupied molecular orbital (HOMO) obtained 
in a self-consistent Hartree-Fock (HF) calculations does 
not account for the observed orientation-dependent ion- 
ization of N2 0| . Only the corresponding length gauge 
MO-SFA 0| and the molecular tunneling theory 0] give 
the correct behavior. 

The orientation dependence of the strong-field ioniza- 
tion of N2 was revisited in Ref. There a simple model 
is made for the 3a g HOMO in terms of a single linear 
combination of atomic orbitals (LCAO) using two atomic 
p orbitals. Then the MO-SFA is applied in the velocity 
gauge to calculate the ionization rate, and it is concluded 
that the velocity MO-SFA, with such an approximation 
for the HOMO, is capable of predicting the observed min- 
imum jjj in the rate for perpendicular orientation of the 
internuclear axis with respect to the linear polarization of 
the laser. Based on these findings it is stated jjj that our 
earlier velocity gauge calculations [2] are inaccurate be- 
cause of inaccuracies in our (i) molecular orbital and (ii) 
generalized Bessel function. Here we refute both these 
claims by giving details of our orbitals and generalized 
Bessel functions. Additionally, we show that apart from 
having the wrong symmetry, the model used for the 3a g 
HOMO of N2 in Ref. [5] is too simple and we therefore 
believe that the agreement in [5] with experiment is ac- 
cidental. To account for the molecular structure more 
LCAOs have to be included as, e.g., in HF calculations. 
With such a self-consistent wave function, the velocity 
gauge MO-SFA, indeed gives the wrong prediction for the 
orientation-dependent ionization Q. Finally, we present 
accurate values from [2j of the generalized Bessel func- 
tion since also in this case some discrepancies exist with 
Ref. (3 - a discrepancy which can be attributed to the 
difference in wavelength (we use 800 nm light Q, in Q 
795 nm light is used 



The idea of the LCAO approach to molecular orbital 
theory is to identify the molecular orbitals in a basis of 
atomic orbitals. For a homonuclear diatomic molecule, 
we write a molecular wave function \£(r) as a superpo- 
sition of atomic orbitals which are centered on each of 
the atomic cores. The nuclei are located at the positions 
±R/2 with respect to the center of the molecule. The 
molecular orbitals are the eigenfunctions of the Fock op- 
erator, i.e., of the kinetic energy operator and the full 
HF potential. These eigenfunctions can be obtained by a 
diagonalization of the matrix representation of the Fock 
operator in the atomic centered basis, (j) n im{r =F R/2). 
A great simplification of the problem is achieved if the 
basis functions are symmetrized according to the molec- 
ular point group (Dqo/i)- Consider, e.g., the atomic Is 
[0ioo (r T orbitals. A change of basis leads to a 

new set of basis functions 

R R 

^a g ls(r) = M ag Is [0100 (*"-—)+ 0100 ( r + y )] U) 

R R 

ipa u is(r) = A/; w i s [0ioo(r - — ) - 0ioo(r + — )]. (2) 

These new basis functions belong to definite symmetry 
classifications in namely E+ ('+' combination) and 

E+ ('-' combination). The 2s, 3s, etc. orbitals are sym- 
metrized similarly. If we extend the basis to contain 
atomic 2p orbitals, the properly symmetrized basis func- 
tions in the E blocks are 

R R 

ip* g 2 P (r) = Afa g 2 P [(t>2io(r - — ) -02io(?" + — )] (3) 

R R 

^a u 2p(r) = A/; w2 p[021o(r - —) + 0210^ + — )], (4) 

for the 2po (%Pz) orbitals, while the and '-' combina- 
tion of the 2p±i are of tt u and 7r g symmetry, respectively, 
and hence not of interest in our present discussion fo- 
cussing on the 3a g HOMO of N2. 

It is important to note that for a quantitative HF calcu- 
lation, the basis is further extended until convergence is 
obtained. We may easily check the correct parity (g or u) 
of the basis function by letting r — > — r and using the fact 
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FIG. 1: (Color online). Planar cuts through the 3a g HOMO 
of N2 obtained in a HF calculation with a (12s7p) /[6s4p] ba- 
sis (see text). Left (right) column: coordinate (momentum) 
space wave function. Upper row: contribution from s type 
basis functions, middle row: contribution from p type basis 
functions, lower row: total HF result in the (12s7p)/[6s4p\ 
basis. 



that the parity of (j) n im is (— We emphasize this point 
since in Ref. [5], the combination of the 2p z orbitals 
is erroneously referred to as the gerade state. Although 
not of relevance for our fully numerical HF calculations, 
we note that in the hydrogenic basis used in Ref. jjj the 
Fourier transform of the antisymmetric combination of 
atomic 2p z orbitals is 



4>* g 2 P (p) = ATa g 2 P sin (p . R/2) 



(1 + V) 3 



Yio(p), (5) 



with Mcr g 2p the normalization constant. We may ver- 
ify the gerade parity symmetry by letting p — > — p and 
using Yi m (— p) — (— l) z F/ m (p). Equation (JJJ is the cor- 
rectly symmetrized version of Eq. (20) of Ref. The 
difference lies in the sine-factor of Eq. © compared with 
a cosine-factor in Ref. [5]. A calculation with the p-type 
basis functions alone would, because of this sine-factor, 
predict suppressed ionization of N2 compared with the 
companion Ar atom - in contrast to experimental and 
theoretical findings Q, 111- Thus having pointed out 
the importance of working with properly symmetrized 
basis function, we now confront the ad hoc single LCAO 
approch Q with our HF calculations. Clearly, but con- 
trary to the approach in Ref. [5], one cannot expect the 



basis functions themselves to be eigenfunctions, instead 
the eigenfunctions are linear combinations of basis func- 
tions. 



III. HARTREE-FOCK CALCULATIONS 



In Ref. we used the GAMESS 0| program for our 
quantitative HF calculations of N2 with a (Ils6p)/[5s3p\ 
Cartesian Gaussian basis The way we specify 

our basis is standard in the quantum chemistry liter- 
ature [nj], but is recalled here for completeness. The 
wave function is expanded in 11 s and 6 p orbitals cen- 
tered on each atom. These basis orbitals all have dif- 
ferent exponentials in the radial parts. The orbitals of 
our HF calculation are obtained by a variation of the ex- 
pansion coefficients until the energy is minimized. This 
procedure is equivalent to an iterative diagonalization of 
matrix representations of the Fock operator. With lit- 
tle loss of accuracy but significant gain in computational 
speed, the optimization procedure can be simplified with 
a contracted basis set. This means that all the expan- 
sion coefficients are not varied independently, contrary 
a fixed relationship is kept between some of the coeffi- 
cients. Such a simplification is used in most quantum 
chemistry software. Five coefficients corresponding to s 
orbitals and three coefficients corresponding to the p or- 
bitals are optimized independently. In order to describe 
the asymptotic behavior we augment the basis by adding 
an extra diffuse s and p orbital. In total, the wave func- 
tion is effectively expanded in a basis of six s and four p 
orbitals on each atom denoted by (12s7p) /[6s4p] basis. 

In Fig. we show results of our calculations for the 
3a g HOMO of N 2 [Fig. [I] (c)] and in its decomposition 
of s [Fig. [I] (a)] and p [Fig. [I] (b)] type basis functions. 
We see that the orbital is composed of both types of basis 
functions. At first sight the p contribution alone seems to 
resemble the orbital quite well. However, in the proxim- 
ity of the nuclei, it is the s contribution that dominates, 
since the p functions have nodes on the nuclei. It is the s 
type contribution that dominates the momentum space 
wave function at low momenta. This fact becomes clear 
in Figs. [I] (d) -(f), where we show the Fourier transforms 
of the functions in Figs.[l](a)-(c). It is clear that both the 
s and p type of basis functions contribute significantly to 
the momentum wave function. Therefore, as discussed 
in Sec. |H] it is insufficient to consider each type exclu- 
sively: The s and p z basis orbitals are not separately 
eigenfunctions of the Fock operator, contrary, the mix- 
ture turns out to be essential. At this point the authors 
of Ref. p| introduce an oversimplification by considering 
the HOMO to be constructed only from a single linear 
combination of two p-orbitals. 
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IV. IMPLICATIONS TO THE MOLECULAR 
STRONG-FIELD APPROXIMATION 

In the MO-SFA one calculates the transition rate from 
the HOMO to a continuum Volkov state through absorp- 
tion of n photons. For a linearly polarized, monochro- 
matic field of frequency u and field strength F , the ve- 
locity gauge result for the angular differential ionization 
rate in atomic units is |l| 



dW n Ar 
— - 2nN e p n 
dp 



*(Pn) (U p -nu) 2 



X J r . 



\F0-P n \ U p 



2u 



(6) 



with N e the number of electrons occupying the HOMO, 
p n — ^2{nuo — U p — I p ) the momentum, U p = Fq/(4cj 2 ) 
the ponderomotive energy, I p the ionization potential of 
the HOMO (0.5725 in the case of N 2 (H), #(p) the mo- 
mentum space wave function and J n (u,v) a generalized 
Bessel function (T^ . 

According to the discussion in Sec. HIH it is important 
to apply the correct momentum space wave function in 
the evaluation of Eq. ©. To illustrate this fact, we con- 
sider the total ionization rate integrated over all angles 
of the outgoing electron and summed over all numbers 
of accessible photon absorptions. We present this to- 
tal ionization rate for different angles between the in- 
ternuclear axis and the polarization axis for the HOMO 
3a g orbital of N 2 . We consider three descriptions of the 
molecular orbital: (i) using s-type basis functions only, 
(ii) using p-type basis functions only and (hi) using the 
self-consistent HF solution. Clearly the approaches (i) 
and (ii) only give a very poor description of the HOMO 
and therefore these results deviate from the result ob- 
tained by using the true HOMO. In Fig. [21 we present 
the results for the three cases at (a) 800 nm Q and (b) 
795 nm [6]. Additionally, we show in Fig.[21the results ob- 
tained with an initial HOMO derived from a grid-based 
HF calculation ^4|. The excellent agreement with the 
(12 si 'p)/ '[6 s4p] basis-state calculation proves the conver- 
gence of the latter approach. We note that panels (a) and 
(b) are very similar and hence the differences between the 
results in Q and Q are not due to the slight difference 
in wavelength. 

From Fig.G](d) we see that the momentum wave func- 
tion in case (i) is nearly spherically symmetric, and hence, 
in Fig. [21 the ionization rate in this case is nearly inde- 
pendent of the molecular orientation. In case (ii) we use 
the momentum function from Fig. ffl(e). Although ex- 
ceptions occur, as evident from Fig.|2{a) below, the gen- 
eralized Bessel functions are typically maximized when p 
is nearly parallel to the polarization direction [Fig.OJb)]. 
In a qualitative analysis, one may therefore expect that 
the maximum rate is obtained, when the polarization axis 
coincides with a direction where |\I/(p)| 2 is large. Corre- 
spondingly, in case (ii) we see from Fig. that the rate 
is maximized when = 0, i.e. when the polarization is 
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FIG. 2: Ionization rate at 2 x 10 14 Wcm 2 of the HOMO of N 2 
(3cr g ) as a function of the angle between the internuclear axis 
and the polarization axis at (a) 800 nm [2] and (b) 795 nm [fj. 
The short-dashed and solid curve are the calculations with the 
HF wave function derived from a grid calculation and with 
the (12s7p)/[6s4p\ basis set in Fig. Q](f), respectively. The 
long-dashed curve with a maximum and the chained curve are 
obtained with the s [Fig. QJd)] and p [Fig. QJe)] contributions 
only, respectively. All rates are normalized to parallel (6 = 0) 
geometry. 



horizontal in Fig.Q(e) and overlaps with a large value of 
the momentum wave function. Finally, we turn to the full 
wave function in case (iii) which is the coherent sum of s 
and p contributions, Fig.[J](f ) . It is the low momenta that 
contribute mostly to the total rate. At these small mo- 
menta, we see from Fig.[J](f) that the absolute square of 
the momentum wave function is maximized along the di- 
rection perpendicular to the molecular axis and therefore 
the rate is maximized with the molecule aligned perpen- 
dicularly to the laser polarization, 6 = it/ 2. This result 
is also shown in Fig. verifying the result of Ref. 0|. 
Note that the rate using the total wave function is not 
the sum of rates of the s and p contributions since their 
contributions to the wave function should be added co- 
herently and with the correct self-consistent amplitudes 
before entering in Eq. We note that the ad- 

mixture of atomic s and p orbitals for the description of 
the HOMO of N2 was also recently pointed out in related 
work on high-harmonic generation [15]. 

In Refs. 0, the values of the generalized Bessel 
function in Eq. © do also not agree. Figure OJa) 
shows our previous result for the square of the general- 
ized Bessel function corresponding to 18 photon absorp- 
tion at 800 nm and at an intensity of 2 x 10 14 W/cm 2 . 
Panel (b) shows the corresponding result at 795 nm jfj. 
Hence, Panel (b) corresponds to Fig. 2(a) of Ref. jjj (since 
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FIG. 3: Polar plot of the generalized Bessel function of order 
n = 18 squared for 2 x 10 14 W/cm 2 , and (a) 800 nm 0, 
(b) 795 nm 0|. The polar angle 
horizontal line. 
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is measured from the 



J- n (u,v) = (—l) n J n (u, —v) pjj) and by comparison we 
see that the results agree. We may therefore attribute 
the difference between the generalized Bessel functions 
in and jjj to stem from the difference in the wave- 
lengths used. It is clear from Fig. that the values of 
the generalized Bessel function are very sensitive to the 
exact value of the arguments. We reproduce the values 
for the generalized Bessel function given in Ref. [16J with 
the input parameters used therein, and we therefore be- 
lieve that our previous and present calculations of the 



generalized Bessel function are correct. Due to the ex- 
cellent agreement at 795 nm with the Bessel functions of 
Ref. [5] we likewise believe that the latter are evaluated 
correctly. 



V. CONCLUSION 

In conclusion, the velocity gauge MO-SFA when based 
on a transition to the continuum from an initial HOMO, 
determined by a HF calculation using standard quantum 
chemistry software, gives the wrong prediction for the 
orientation-dependent rate. In this approximation only 
the length gauge MO-SFA (and the molecular tunneling 
theory) [2] give the observed minimum [3] in the rate for 
the perpendicular geometry. This conclusion contrasts 
the findings in a recent work p]. In that work, however, 
only a single LCAO of two p-orbitals (of wrong symme- 
try) was considered for the description of the HOMO and 
no self-consistent wave function was applied. From Fig. [21 
we clearly see, that an application of a too simple wave 
function, e.g., one with only p basis functions, leads to 
a prediction accidentally in qualitative agreement with 
experiment [3j. If one nevertheless would be encouraged 
by this agreement to believe that the single LCAO of p 
orbitals describes the physics well, one should note that 
such a wave function due to the presence of the sine- 
factor in (JU would predict suppressed ionization in N2 
compared with Ar which is in contrast with experiments 
and theory 0, 0, . In closing we note that also in the 
atomic case, the length gauge version of the SFA was re- 
cently shown to be superior to the velocity gauge version 
when compared with results obtained by integrating the 
time-dependent Schrodinger equation [l7j . 
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